Computing flow past a cylinder in ViscousFlow


In [1]:
using ViscousFlow

In [2]:
using Plots
pyplot()
clibrary(:colorbrewer)
default(grid = false)

Solve flow past a stationary body

Set the flow parameters


In [3]:
Re = 200; # Reynolds number
U = 1.0; # Free stream velocity
U∞ = (U,0.0);

Set up points on the body.


In [4]:
n = 100;
body = Bodies.Ellipse(0.5,n)


Out[4]:
Elliptical body with 100 points and semi-axes (0.5,0.5)
   Current position: (0.0,0.0)
   Current angle (rad): 0.0

Transform the body to lie at a specific position. Here, we just put it at the origin.


In [5]:
cent = (0.0,0.0)
α = 0.0
T = RigidTransform(cent,α)
T(body) # transform the body to the current configuration


Out[5]:
Elliptical body with 100 points and semi-axes (0.5,0.5)
   Current position: (0.0,0.0)
   Current angle (rad): 0.0

Set up the domain and plot the body inside it.


In [6]:
xlim = (-1.0,3.0)
ylim = (-1.0,1.0);

In [7]:
plot(body,xlim=xlim,ylim=ylim)


Out[7]:

Now set up the coordinate data for operator construction


In [8]:
X = VectorData(body.x,body.y);
#X̃ = VectorData(body.x̃,body.ỹ);

Set the grid cell size and time step size.


In [9]:
Δx = 0.02
cfl = 0.5
Δt = min(cfl*Δx,0.5*Δx^2*Re)


Out[9]:
0.01

Now set up the system

Set up the state vector and constraint force vector for a static body


In [10]:
sys = Systems.NavierStokes(Re,Δx,xlim,ylim,Δt,U∞ = U∞,  = X, isstore = true)


Out[10]:
Navier-Stokes system on a grid of size 202 x 102

In [11]:
w₀ = Nodes(Dual,size(sys));
f = VectorData(X);

In [12]:
wf = Systems.PointForce(w₀,(1.5,0.0),10.0,4.0,1.0,sys)


Out[12]:
Transient point force applied on the Nodes{Dual,202,102} field.
   strength = 10.0
   location = (1.5, 0.0)
   central time = 4.0
   half-interval = 1.0

In [13]:
xg, yg = coordinates(w₀,dx=Δx,I0=Systems.origin(sys))


Out[13]:
(-1.01:0.02:3.0100000000000002, -1.01:0.02:1.01)

Set up the integrator here


In [14]:
plan_intfact(t,u) = Systems.plan_intfact(t,u,sys)
plan_constraints(u,t) = TimeMarching.plan_constraints(u,t,sys)
r₁(u,t) = TimeMarching.r₁(u,t,sys) + wf(t)
r₂(u,t) = TimeMarching.r₂(u,t,sys)

@time ifherk = IFHERK(w₀,f,sys.Δt,plan_intfact,plan_constraints,(r₁,r₂),rk=TimeMarching.RK31,isstored=true)


  6.683975 seconds (12.66 M allocations: 1.085 GiB, 6.29% gc time)
Out[14]:
Order-3 IF-HERK integrator with
   State of type Nodes{Dual,202,102}
   Force of type VectorData{100}
   Time step size 0.01

Initialize the state vector and the history vectors


In [15]:
t = 0.0
u = deepcopy(w₀)

fx = Float64[];
fy = Float64[];
thist = Float64[];

uhist = [];
tsample = 0.2; # rate at which to store field data

Advance the system!

Set the time range to integrate over.


In [16]:
tf = 30.0;
T = Δt:Δt:tf;

In [17]:
@time for ti in T
    global t, u, f = ifherk(t,u)
    
    push!(thist,t)
    push!(fx,sum(f.u)*Δx^2)
    push!(fy,sum(f.v)*Δx^2)
    (isapprox(mod(t,tsample),0,atol=1e-8) || isapprox(mod(t,tsample),tsample,atol=1e-8)) ? push!(uhist,u) : nothing
end
println("solution completed through time t = ",t)


254.074707 seconds (1.12 G allocations: 54.627 GiB, 4.88% gc time)
solution completed through time t = 30.00000000000189

In [18]:
thist[end]


Out[18]:
30.00000000000189

Plotting

Basic plot


In [19]:
plot(xg,yg,uhist[end],levels=range(-0.1,stop=0.1,length=30), color = :RdBu,xlim=(-1+Δx,3-Δx),ylim=(-1+Δx,1-Δx),width=1)
plot!(body)


Out[19]:

Make a movie


In [20]:
anim = @animate for i = 1:length(uhist)
    plot(xg,yg,uhist[i],levels=range(-0.1,stop=0.1,length=30),color = :RdBu,xlim=(-1+Δx,3-Δx),ylim=(-1+Δx,1-Δx))
    plot!(body)
end
gif(anim,"cylinderRe200.gif")


┌ Info: Saved animation to 
│   fn = /Users/jeff/ViscousFlow/examples/cylinderRe200.gif
└ @ Plots /Users/jeff/.julia/packages/Plots/qh1wV/src/animation.jl:90
Out[20]:

In [21]:
plt = plot(layout = (2,1), size = (600, 400))
plot!(plt[1],thist,2*fy,xlim=(0,30),ylim=(-2,2),xlabel="Convective time",ylabel="\$C_L\$",legend=false)
plot!(plt[2],thist,2*fx,xlim=(0,30),ylim=(0,4),xlabel="Convective time",ylabel="\$C_D\$",legend=false)
plt


Out[21]:

In [22]:
using Statistics

In [23]:
thist[3000]


Out[23]:
30.00000000000189

In [30]:
Statistics.mean(2*fx)


Out[30]:
1.307759701804912

In [32]:
using Pkg

In [33]:
Pkg.add("FFTW")


  Updating registry at `~/.julia/registries/General`
  Updating git-repo `https://github.com/JuliaRegistries/General.git`
 Resolving package versions...
  Updating `~/.julia/environments/v1.0/Project.toml`
  [7a1cc6ca] + FFTW v0.2.4
  Updating `~/.julia/environments/v1.0/Manifest.toml`
 [no changes]

In [34]:
using FFTW

In [35]:
fft(fx)


Out[35]:
3000-element Array{Complex{Float64},1}:
  1961.6395527073682 + 0.0im               
  -130.9955220894658 + 154.9277919290999im 
  -92.35257168263149 + 30.194809081364568im
 -124.98358692223584 + 4.878021444668im    
 -130.52957799712146 + 7.502131878946741im 
 -136.84704259473472 + 6.267905305803859im 
 -135.19963408568682 + 8.872988300622477im 
  -135.9271764065013 + 7.7679493768210754im
 -136.51605357994052 + 7.779823718421892im 
 -138.73147306418542 + 11.674193041482201im
 -128.62169828772096 + 13.861622725189225im
 -132.64901382568348 + 1.77211738089636im  
   -135.732155895115 + 8.242809926055713im 
                     ⋮                     
   -135.732155895115 - 8.242809926055713im 
 -132.64901382568348 - 1.7721173808963535im
   -128.621698287721 - 13.861622725189221im
  -138.7314730641854 - 11.674193041482201im
  -136.5160535799405 - 7.779823718421893im 
 -135.92717640650127 - 7.767949376821077im 
  -135.1996340856868 - 8.872988300622477im 
 -136.84704259473475 - 6.267905305803859im 
 -130.52957799712146 - 7.502131878946741im 
 -124.98358692223584 - 4.8780214446680015im
  -92.35257168263153 - 30.194809081364575im
 -130.99552208946582 - 154.9277919290999im 

In [36]:
nt = length(thist)
nfreq = (nt-1)÷2+1
fsamp = 1/Δt
freq = 0.5*fsamp*range(0,stop=1.0,length=nfreq)
fhat = fft(fy[1:nt-1])/(nt-1)
plot(freq,4*abs.(fhat[1:nfreq]),xlim=(0,1))


Out[36]:

In [39]:
xg,yg = coordinates(u,dx=Δx,I0=Systems.origin(sys))
ψ = deepcopy(uhist[end])
ψ .= -(sys.L\uhist[end])*Δx .+ repeat(collect(yg)',size(sys,1))
plot(xg,yg,ψ,trim=1,clim=(-0.1,0.1))
plot!(body,fillcolor=:black,linecolor=:black)


Out[39]:

In [ ]: